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ABSTRACT 

Asteroseismology is able to conduct studies on the interiors of solar-type stars from 
the analysis of stellar acoustic spectra. However, such an analysis process often has to 
rely upon subjective choices made throughout. A recurring problem is to determine 
whether a signal in the acoustic spectrum originates from a radial or a dipolar oscil- 
lation mode. In order to overcome this problem, we present a procedure for modelling 
and fitting the autocovariance of the power spectrum which can be used to obtain 
global seismic parameters of solar-type stars, doing so in an automated fashion with- 
out the need to make subjective choices. From the set of retrievable global seismic 
parameters we emphasize the mean small frequency separation and, depending on the 
intrinsic characteristics of the power spectrum, the mean rotational frequency split- 
ting. Since this procedure is automated, it can serve as a useful tool in the analysis of 
the more than one thousand solar-type stars expected to be observed as part of the 
Kepler Asteroseismic Investigation (KAI). We apply the aforementioned procedure 
to simulations of the Sun. Assuming different apparent magnitudes, we address the 
issues of how accurately and how precisely we can retrieve the several global seismic 
parameters were the Sun to be observed as part of the KAI. 

Key words: methods: data analysis - methods: statistical - stars: oscillations. 



INTRODUCTION 

Seismology of solar-type oscillators is a powerful tool that 
can be used to increase our understanding of stellar struc- 
ture and evolution. Oscillations in main-sequence stars and 
subgiants have been measured thanks to data collected from 
grou nd-based high-precision sp ectroscopy (for a review see 
e.g.. iBedding fc Kieldsenll200a ) and, more recently, to pho- 
tometric space-based mi s sions such as CoR oT (see e.g., 
lAppourchaux et alj [20ol : iMichel et al.l l200&t ). The Kepler 
mission (for a discussion on the expected results of the 
asteroseismic investigat ion see IChristensen-D alsgaar d et al.l 
l200Sl : lKaroff et al.ll2009T ) will lead to a revolution in the field 
of asteroseismology of solar-type oscillators, since it will in- 
crease by more than two orders of magnitude the number 
of stars for which high-quality observations will be avail- 
able, while allowing for long-term follow-ups of a selection 
of these targets. As of the time of writing of this article, first 
results arising from the Kepler aste roseismic programme 
had already been made available (B edding et al.l l20ld : 



Chapli n et al. |2010|: iGilliland et alj|20ld: Grigahcene et al 
201C : lHekker et aljhoiCtlKolenberg et al.ll2oTd ; IStello et a! 
201O) 



Due to the large number of stars observed with Ke- 
pler, automated and innovative analysis pipelines/tools 
are needed in order to cope with the plenitude of avail- 



able data (see e.g.. iHekker et al. 


120091: iHuber et al.l 


2009 


Kallinger et al.l 20101: Karoff et al. 


20ld; Mathur et all 


201C 


Mosser & Aouourchaux 2009: Roxburgh 2009). The 


auto- 



mated pipelines dedicated to the analysis of acoustic spec- 
tra that have been developed so far aim mainly at measuring 
the frequency of maximum amplitude, v max , the maximum 
mode amplitude, A max , and the mean large frequency sep- 
aration, Ai/. In the present study we give continuity to this 
work by presenting a tool capable of modelling and fitting 
the ACP^j (AutoCovariance of the Power Spectrum) of a 
solar-type oscillator. The current version of the tool accepts 
as free model parameters the mean small frequency sepa- 
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1 To be precise, we compute the autocovariance of the power 
density spectrum. 
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ration parameter, Do (see £13.1.21 for the definition of this 
parameter), the mean rotational frequency splitting, i/ B , the 
mean linewidth, (F), and the stellar inclination angle, i. The 
output generated by this tool will hopefully contribute to 
further explore the diagnostic potential of solar-like oscil- 
lations, e specially that of the large a nd small separations 
(see e.g., IChristensen-Dalsgaardl [20041 ). These two quanti- 
ties can in turn be used to provide estimates of the radii, 
mas ses and ages of sola r- like stars with consistent uncertain- 
ties (|Karoff et al.ll2010l ). 

We start in [|2] by providing a logical basis that led us 
to the development of this particular analysis tool. This is 
followed by a thorough description of its implementation in 
Sj3] In Sj4]we apply the tool to simulated data, more specif- 
ically, to a solar analog as it would be seen by Kepler. A 
discussion and conclusions are presented in EJS] 



2 RATIONALE BEHIND THE 

ACPS-MODELLING PROCEDURE 

Unlike a canonical peak-bagging p rocedure (see e.g., 
I Anderson. Duvall Jr. fc Jefferieslll990l ). in the framework of 
the ACPS-modelling procedure individual modes do not 
need to be tagged by angular degree, i.e. by wave number 
I. This is required in peak-bagging to ensure the correct 
model is fitted to the observed modes. Instead we build a 
model that describes both the global and average properties 
of the most prominent modes, which are in turn encoded on 
a small set of free model parameters. This might be regarded 
as a major advantage, especially if we recall the ambiguity 
found in mode (angular degree) identification concerning the 
CoRoT target HD 49933, only recentl y solved after a new 
longer time series was ma de available |Appourchaux et al.l 
|200S| ; iBenomar et~ai1l2009l ) . 

As it is currently implemented, the ACPS-modelling 
procedure assumes the small frequency separation to be 
constant and not a function of frequency. Especially for 
evolved stars this was thought n ot to be a good approxi- 
mation (|Soriano fc Vauclairll2008i ). H owever, in the light o f 
the first Kepler results on red giants (|Bedding et al.l [2010^. 
it is perfectly valid to assume - at least for low-luminosity 
red giants - a mean value of the small frequency separa- 
tion between adjacent modes with £ = and 1 = 2, 5vo2, 
it being in fact a nearly constant fraction of Ais. There- 
fore, the procedure presented herein can in principle also 
be employed in the case of evolved stars displaying solar- 
like oscillations. The ACPS-modelling procedure will how- 
ever not provide sensible output when m ixed modes (see 
e.g., lAizenman. Smevers fc WeigertJ Il977h are present. In 
any event, the presence of mixed modes will manifest clearly 
in the ACPS, either as a significant broadening of the peaks 
or by the appearance of extra peaks. 

The way information is presented in the ACPS makes 
it very amenable for fitting. In fact, the ACPS of a solar- 
type oscillator will show prominent features at multiples of 
half the large frequency separation, a clear manifestation of 
the regular frequency structure of the acoustic spectrum. 
The shape of these features will depend upon the Svo2 and 
Si/oi spacings (see H3.1.2l for a definition of the latter), mode 
lifetime, rotation and stellar inclination. Fig. [1] displays an 
artificial acoustic spectrum of a solar analog (top panel) as 
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Figure 1. Top Panel: Artificial acoustic spectrum of a solar ana- 
log evidencing a regular frequency structure. Bottom Panel: Fit 
(black solid line) to the corresponding ACPS obtained using the 
MAP as our choice of summary statistic to represent the param- 
eter posterior distributions. In order to enhance the visibility of 
the features in the ACPS we have slightly smoothed it by run- 
ning a boxcar average (this is done both to the observed and 
to the model ACPS). Notice that the artificial star (Boris) has 
Af~135 /iHz (see discussion in Q. 

it would be obtained from a 30-day long Kepler's time series 
together with a fit to its ACPS (bottom panel). Notice the 
prominent features in the ACPS around 135 /iHz and half 
that value. The latter feature is clearly split into two peaks 
which constitutes a signature of Suqi . A signature of the 5vq2 
spacing is also present in the wings of the feature at Av. The 
widths of either feature are generally a signature of the mode 
lifetime of the oscillations, whereas the presence of any fine 
structure is a signature of rotation and stellar inclination. 
The ACPS depicted in Fig. [T] does indeed provide us with 
information about the large and small frequency separations, 
and mode lifetime. The question we address first is thus: 
How do we extract this information from the autocovariance 
of the power spectrum ? 



3 THE ACPS-MODELLING PROCEDURE 

We intend to fit the observed ACPS to an ACPS created 
from a model of the p-mode power density spectrum (PDS) 
that is described by a slightly modified asymptotic relation, 
which in turn depends only on a few free parameters (JO]). 
Although we do not fully explore its potential, we opt for a 
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Bayesian approach to the fitting problem for the simple rea- 
son that any prior information might be taken into account, 
particularly useful when constraining the parameter space 
( . Most importantly, we aim at obtaining a full joint 
posterior probability density function (PDF) concerning the 
set of model parameters and hence we end up employing a 
Markov chain Monte Carlo (MCMC) sampler (S£OJ. 



3.1 Modelling the observed ACPS 

The ACPS is computed as a function of the frequency lag, 
L, according to the following expression: 



ACPS(L) = 



N 



E 

3=0 



(Pj - P) {P ]+ L - P) , 



(1) 



where Pj is the particular value taken by the PDS at the 
fixed frequency bin j, P is the mean value of the sample 
population {Pj}, i.e. the mean value of the PDS within the 
frequency interval of interest, and N is the total number of 
bins within this same frequency interval. The frequency in- 
terval of interest should be in accordance with the p-mode 
range. The lag, L, which is expressed in terms of number of 
frequency bins, is allowed to vary between the corresponding 
frequency values of 10 /iHz and 4/3 Av. The reason why we 
constrain L to this range has to do with the fact that we are 
primarily interested in measuring the small frequency sepa- 
ration. The signature of the small frequency separation is in 
fact mainly conveyed by the features at Av/2 and Av. For 
L less than 10 fiRz the ACPS is dominated by the signature 
of the window function, whereas for L greater than approx- 
imately 4/3 Av it reproduces the features seen below this 
value of the lag and thus becoming redundant. Finally, the 
ACPS is normalised after division by its maximum value. 

As already stated, modelling the observed ACPS con- 
sists of two steps: (i) building a model of the subjacent 
p-mode PDS and (ii) computing its autocovariance as ex- 
plained above. In what follows we provide the reader with 
the basic theory behind the power spectrum of solar-like 
oscillations and end up describing how the model p-mode 
spectrum is generated. 



3.1.1 The model 

Solar-like oscillations are global standing acoustic waves, 
also known as p m odes, driven by n ear-surface turbulent 
convection (see e.g.. lBalmforthiri992al lbT). Radial oscillation 
modes are characterised by the radial order n, whereas non- 
radial modes are additionally characterised by the non-radial 
wave numbers I and m. 

Ignoring any departure from spherical symmetry, non- 
radial modes differing only on the azimuthal wave number m 
are degenerate. Stellar rotation removes the (2£+ l)-fold de- 
generacy of the frequency of oscillation of non-radial modes. 
When the angular velocity of the star, Q, is small and in 
the case of rigid-body rotation, the frequency of a (n, i, m) 



mode is given to first order by (|Ledou: 



3qucncy 
iHII): 



(2) 



Vnlm = Vnl + J71— (1 — C n l) . 

Ztv 

The kinematic splitting, mfi/(27r), is corrected for the ef- 
fect of the Coriolis force through the dimensionless quan- 



tity Cni > 0. In the asymptotic regime, i.e. for high-order, 
low-degree p modes, rotational splitting is dominated by ad- 
vection and the splitting between adjacent modes within a 
multiplet is u B ~ Sl/(2ir). Typically, v s q =0.4 fiRz. 

The PDS of a single mode of oscillation is distributed 
around a li mit spectrum with an exponential prob ability 
distribution l|Woodard|[l98i : iDuvall fc Harvey^ ! 986). This 
limit spectrum contains the information on the physics of 
the mode and may be described as a standard Lorentzian 
profile near the resonance. It follows that the overall lim it 
p-mode spectrum is given by (see e.g.. iFletcher et alj|2006h : 



Vnl j l^s i J- n£m j 



n=riQ £—0 m — — £ 1 -f- 



2 (v— v n i — mv s ) 



(3) 



where Sntm is the mode height, v n t is the central frequency 
of the (n, £) multiplet and T n t m is the mode linewidth. In 
the case of solar-type oscillators and for low angular degree 
£, we can assume that T is a function of frequency alone. V is 
related to the mode lifetime, r, through r = (7rr) _ . JV„ de- 
scribes the background signal originating from both granula- 
tion and activity. Notice that we are assuming that a mode is 
uncorrelated with any other modes or with the background 
signal. The stellar background signal is commonly modelled 
as a sum of power laws describing these physical phenomena 
(|Harvevlll985l : lAigrain. Favata fc Gilmorell2004h : 



N, 



M = E T 



+ (2nB k v) c 



+ N, 



(4) 



{Ak} and {P>k} being, respectively, the corresponding am- 
plitudes and characteristic time scales, whereas the {Ck} are 
the slopes of each of the individual power laws. A flat offset, 
N, is needed in order to model photon shot noise. 



3.1.2 The mode frequencies 

Before we start modelling the observed ACPS we first need 
to select a suitable frequency interval for the purpose of 
our study. This interval should coincide with the frequency 
range in which p modes are located - i.e. from the frequency 
of the fundamental (n=0) mode, vj, up to the atmospheric 
acoustic cut-off frequency, i^ ac . 

We compute an estimate of the frequency of the fun- 
damental mode by assuming that it scales with the mean 
stellar density and thus with the mean large frequency sep- 
aration: 



v { = v m (Av/Avq) . 



(5) 



where V{q = 258 ^iHz and A^q — 135 ^iHz. 

An estimate of the atmospheric acoustic cut-off fre- 
quency is calculated by assumin g that it scales with the fre - 
quency of maximum amplitude dKieldsen fc BeddingHl995T ): 



(^max / ^maxy ) ■ 



(6) 



For the Sun one has f aC Q = 5300 fiRz and ^ ma x© = 3100 piHz. 

We proceed with the assignment of frequency values to 
modes with degree £ = 0, 1, 2 (see H3.1.3l for an explanation of 
this upper lim it on £) according to the asymptotic relation 
l|Tassoullll980l ): 
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Vni = Av (n + -I + e) - l{l + 1)D , 



(7) 



where the quantity e is sensitive to the surface layers and 
was taken to be 1.5. Moreover, Do is sensitive to the sound 
speed gradient near the core. If Eq. (0 holds exactly then 
it follows tha10 Svo2 = 6Do, Si/%s = 10-Dp and Svoi = 2Dq. 

The automated Kepler pipeline (He kker et al.l I2009T I 
used to supply input for the procedure described herein gen- 
erates as output the smooth second order change in the large 
frequency separation as a function of n, dAu/dn. In cases 
where an estimate of dAu/dn is available we opt for a mod- 
ified version of Eq. (0, which includes a second order cor- 
rection: 

v nl = Au(n+^£ + e)-£(£ + l)D a + (n-n m ^) 2 AAv J An .( 8 ) 

The overtone with the highest power, n max , is given bjf] 
round[(^ max — Vf)/Av\. 

For the sake of clarity we should again stress that indi- 
vidual modes in the observed PDS do not need to be tagged 
by angular degree £ prior to the fit. Assignment of frequency 
values to modes in the model of the PDS, according to either 
Eq. (0 or Eq. (0, is in fact done in complete ignorance of 
the correct mode tagging. 



3.1.3 The mode heights 

The height of a multiplet component can be expressed as: 



Sntm — £tm{i) S n £ — £l m (i) Vl O-nl ■ 



(9) 



The £tm(i) factor represents mode visibility within a multi- 
plet and i is the inclination angle between the direction of 
the stella r rotation axis and the line of sight. This factor is 
given by (jGizon fc Sorankill2003h : 



(£-\m\)\ [ H .J 2 



(10) 



where P™(x) are the associated Legendre functions. The 
quantity V 2 is an estimate of the geometrical visibility of the 
total power in a (n, £) multiplet as a function of £, whereas 
ctn,e depends mainly on the frequency and excitation mech- 
anism. Eq. ((9| is only strictly valid under one assumption: 
When the stellar flux is integrated over the full apparent 
disc, one must assume that the weighting function depends 
only on the distance to the disc centre. In this case, the ap- 
parent mode amplitude can effectively be separated into two 
factors: £i m (i) and V 2 . This assumption holds very well in 
the case of intensity measurements, since the weighting func- 
tion is then mainly linked to the limb-darkening, whereas 
for velocity measurements, departures might be observed 
due to asymmetries in the velocity field induced by rot a- 
tion (see lBallot. Garcia fc LambertllioOr] : iBallot et ai.Jl200St 
and references therein). The heights of non-radial modes are 
commonly defined based on the heights of radial modes ac- 
cording to Eq. j9l), and takin g into account the Ve /Vo ratios 
(see e.g.. lBedding et al]|l996l ). Note that £ = modes consti- 
tute a sensible reference since they are not split by rotation. 



2 <5fi3 is the spacing between adjacent modes with 1=1 and I = 3, 
and Svoi is the amount by which 1 = 1 modes are offset from the 
midpoint between the 1 = modes on either side. 

3 The round[ ] function rounds the argument to its closest integer. 



In the present case of stellar observations or when observing 
the Sun- as- a- star, the associated whole-disc light integra- 
tion strongly suppresses high-degree modes due to the lack 
of spatial resolution. Stellar observations are hence mostly 
sensitive to high-order acoustic eigenmodes with £ ^ 3. 

We compute the power envelo pe for radial m o des a s 
a function of frequency according to iKieldsen et al l l|200Sft . 
We start by subtracting a fit to the background signal from 
the observed PDS. The residual spectrum thus obtained is 
heavily smoothed over the range occupied by the p modes by 
convolving it with a Gaussian having a FWHM of 4 Av. Fi- 
nally, we multiply the smoothed, residual spectrum by Av/c, 
where c is a factor that measures the effective number of 
modes per order and taken to be 3.03. The height, S„o, of 
a radial mode in un its of power density is then given by 
jChaplin et alJl200Sh : 



2A 2 T 

TTTFnO + 2 



(11) 



where A 2 is the total power of the mode (as determined from 
the power envelope for radial modes) and T is the effective 
observational length. 



3.1.4 Setting up the model 

In order to generate the model p-mode spectrum, our mod- 
elling procedure requires the following input: The observed 
PDS sampled at the true resolution (i.e. computed over a 
grid of uncorrelated frequency bins), a fit to the background 
signal, fmax an d Av. The gradient of the large frequency 
separation with n, dAv/dn, is optional. This input, with 
the exception of the observed PDS, might be obtained from 
the automated Kepler pipelines. 

A common linewidth value is assigned to all the modes 
(a sensible assumption if the range in frequency of the spec- 
trum being tested is not too wide). The fitted parameter (r) 
might then be interpreted as a mean linewidth. 

The model p-mode spectrum (including the background 
signal) is finally assembled according to Eq. (J3J), taking into 
account the global seismic parameters fitted in our proce- 
dure, namely, Do, ^s, (F) and i. In case we would want to 
allow for the effect of the windowing and statistical weight- 
ing, then we should at this stage convolve the model p-mode 
spectrum with the spectral window. 



3.2 Bayesian inference 

3.2.1 Parameter estimation 

In the framework of Bayesian parameter estimation, a par- 
ticular model of the observed ACPS is assumed to be true 
and the hypothesis space of interest relates to the values 
taken by the model parameters, = {-Do, v B , (P), i}. These 
parameters are continuous, meaning that the quantity of in- 
terest is a PDF, in contrast to traditional point estimate 
methods. Let us state Bayes' theorem: 



P (&\D,I) 



(&\I)p(D\&,I) 



(12) 



p(D\I) 

where D represents the available data and the prior infor- 
mation is represented by /. p(&\I) is called the prior prob- 
ability, whereas p(&\D, I) is called the posterior probability. 
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The likelihood function is represented by p(D\®, I). p(D\I), 
the global likelihood, simply plays the role of a normalisation 
constant. 

The procedure of marginalisation makes it possible to 
derive the marginal posterior PDF for a subset of parameters 
0a, by integrating out the remaining parameters 0b, the 
so-called nuisance parameters: 

p(® A \D,I) = J p(& A ,& B \D,I)d& B . (13) 

Furthermore, assuming that the prior on ©a is independent 
of the prior on the remaining parameters, then by applying 
the product rule we have: 

p(0a,0b|/)=p(0a|/)p(0b|0a,/)=p(0a|/)p(0b|/).(14) 

The main advantage of the Bayesian framework when 
compared to a frequentist approach is the ability to incorpo- 
rate relevant prior information through B ayes' theorem and 
to evaluate its effect on our analysis (see e.g.. lBrewerll200§l . 
for possible applications of Bayesian probability theory in 
Astrophysics). 



3.2.2 The likelihood function 

We would like to specify the likelihood function, i.e. the 
joint PDF for the data sample. We know that, for a given 
frequency lag Li € [I m i n ,L max ]: 



3.3 Markov chain Monte Carlo (MCMC) 



ACPSo(Li) = ACPS m (Li) + e, 



(15) 



where ACPS and ACPS m are the observed and model 
ACPS, respectively. The error term e, follows a GaussiarQ 
distribution with zero mean and standard deviation o~i. 
Hence, assuming that the model is deterministic, i.e. true, 
we can write 



/[ACPSo(Li)] 



exp 



2cr 2 



exp 



aV2n 

ACPSp^) - ACPS m (L 
2a 2 



(16) 



where / [ACPS (£i)] is the probability density that the ob- 
served ACPS takes a particular value ACPS (ii) at a fre- 
quency lag Li. Notice that we are assuming cr, to be constant 
over the whole range spanned by Li. 

By ignoring the effect of correlation between points in 
the ACPS we arrive at the following expression for the like- 
lihood function: 



£(0) = p(D\®, /) = n/ [ACPSo(Li) 



(17) 



While the effect of correlation clearly must be present, the 
fact that we ignore it can be shown to only affect attempted 
error calculations an d not the fitted values themselves (see 
iFletcher et al.l 12006. and references therein). This does not 
constitute a problem in the present case, since we will be 
employing a MCMC sampler in order to obtain the marginal 
posteriors for each of the model parameters, from which we 
can estimate the uncertainties. 



4 Each point in the ACPS is the sum of a large number of in- 
dependent and identically distributed random variables. There- 
fore, the Central Limit Theorem applies and the distribution of 
ei should approach the normal distribution. 



After inspection of Eq. (|13[) . the need for a mathematical 
tool that is able to efficiently evaluate the multi-dimensional 
integrals required in the computation of the marginal poste- 
riors becomes clear. This constitutes the principal rationale 
behind the method known as Markov chain Monte Carlo, 
which aims at drawing samples from the target distribution, 
p(&\D, I), by constructing a pseudo-random walk in model 
parameter space. Such a pseudo-random walk is achieved by 
generating a Markov chain, whereby a new sample, ©t+i, 
depends solely on the previous sample, ©t, in accordance 
with a time-independent quantity called the transition ker- 
nel, p(® t +i\&t): 



g(0 t+ i|0 t ) min 



p(e t +i|o*) 

p(® t+1 \D,I) q(® t \® t+1 ) 



p(®t\D,I) q(® t+1 \® t ) 



(18) 



where q(®t+i|®t) is a proposal distribution centred on f . 
After a burn-in phase, p(© t +i|©t) is able to generate sam- 
ples of with a probability density converging on the 
target distribution. The algorithm that we employ in or- 
der to generate a Mark ov chain was initially proposed by 
Metropolis et alj (|l953T l. and subsequently generalised by 
Hastings! (|l970h . this latter version being commonly referred 
to as the Metropolis-Hastings algorithm. 

The Metropolis-Hastings MCMC algorithm just out- 
lined above runs the risk of becoming stuck in a local maxi- 
mum of the target distribution, thus failing to fully explore 
all regions in parameter space containing significant prob- 
ability. A way of overcomi ng this difficulty i s to employ 
parallel tempering (see e.g.. lEarl fc Deerr]|2005r i. whereby a 
discrete set of progressively flatter versions of the target dis- 
tribution is created by introducing a temperature parameter, 
T . In practice, use is made of its reciprocal, /3 = 1/T, referred 
to as the tempering parameter. By modifying Eq. (|12p . we 
generate the tempered distributions as follows: 

P (G\D,p,I) =Cp{®\I)p{D\@,lf , < /3 ^ 1, (19) 

where C is a constant. For /3 = 1, we retrieve the target 
distribution, also called the cold sampler, whereas for /3<1, 
the hotter distributions are effectively flatter versions of the 
target distribution. By running such a set of cooperative 
chains in parallel and by further allowing for the swap of 
their respective parameter states, we enable the algorithm to 
sample the target distribution in a way that allows for both 
the investigation of its overall features (low-/? chains) and 
the examination of the fine details of a local maximum (high- 
/3 chains). See Appendix [Al for the current version of the 
parallel tempering Metropolis-Hastings algorithm written 
in pseudocode. 

Moreover, based on a stati stical control sy stem (CS) 
similar to the one described in ICregorvl (|20051), we auto- 
mate the process of calibration of the Gaussian proposal 
crcs-values, which specify the direction and step size in pa- 
rameter space when proposing a new sample to be drawn. 
The optimal choice of {ercs} is closely related to the aver- 
age rate at which proposed state changes are accepted, the 
so-called acceptance rate. The control system makes use of 
an error signal to steer the selection of the crcs-values dur- 
ing the burn-in stage of a single parallel tempering MCMC 
run, acting independently on each of the tempered chains. 
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The error signal is proportional to the difference between 
the current acceptance rate and the target acceptance rate. 
As soon as the error signal for each of the tempered chains 
is less than a measure of the statistical fluctuation expected 
for a zero mean error, the control system is turned off and 
the algorithm switches to the standard parallel tempering 
MCMC. 



4 APPLICATION TO SIMULATED DATA: THE 
SOLAR TWIN BORIS 

We display the results obtained when applying the ACPS- 
modelling procedure to the particular case of Boris, an ar- 
tificial main-sequence star created in the framework of the 
AsteroFLAG group for the purpose of conducting hare and 
hounds. Briefly, Boris is what we might call a solar analog 
(L = 1.0 L w ; T eff = 5778 K; R = 1.00 R w ). W e refer the 
reader to IChaplin et all (|2008h and IStello et all (|2009h for 
further insight into AsteroFLAG's activities. 




4.1 Experimental setup 

Boris was assumed to have an apparent visual magnitude 
of either V = 8 or V = 9 or V — 10 (these values corre- 
spond to the bright end of the nominal apparent magnitude 
target range for Kepler). In each case, 250 realisations of a 
PDS were generated. The goal was to assess both the ac- 
curacy and the precision of the several Bayesian summary 
statistics of the fitted parameters, namely, the MAP (Maxi- 
mum A Posteriori) , the median, the mean and the marginal 
posterior mode. These power density spectra were directly 
generated in Fourier space and include contributions aris- 
ing from p modes, granulation, activity, photon shot noise 
and instrumental noise. They correspond to 30-day long 
time series with a 60-second cadence, as it would be ex- 
pected for Kepler's high-cadence survey targets. The stellar 
inclination angle, i, entering the simulations was allowed to 
vary and was drawn from a probability distribution that 
assumes random orientations - sin(i). Reference values are 
given for the other model parameters by: Dfj ef = 1.43 /iHz, 
vl cl = ^ 80 = 0.4 lifiz and (r) rcf = 1.95 rtHz. All the remaining 
relevant quantities that enter the simulations were assigned 
'solar' values. 

A prototype of the Birmingham -Sheffield Hallam au- 
tomated pipeline |Hekker et al.l 20091 ') was then run on the 
simulated power density spectra in order to retrieve useful 
input for the ACPS-modelling procedure, to be specific, a fit 
to the background signal, Av and dA^/dn. The frequency 
interval of interest (for purposes of the computation of the 
ACPS) ranges from 2100 to 4500 /iHz and f ma x has been 
fixed at 3300 /iHz. 

Figs. [T] [2] and [3] display the typical graphical output of 
the ACPS-modelling procedure resulting from the analysis 
of a single realisation of the PDS of Boris (with V = 8 and 
a random orientation characterised by i = 32.4°). Note that 
throughout the analysis of the full set of realisations of the 
power spectrum, uniform priors have been imposed on Do, 
v s and (r), whereas p[i\I)=sm.{i). 



Figure 2. Sequence of 1 X 10 5 samples (at a thinning interval of 
1) drawn from the target distribution (/3 = 1) by a Metropolis- 
Hastings Markov chain Monte Carlo. Top Panels: Behaviour 
of log 10 [p(0|7) £(©)] as a function of the iteration number (a 
means of visually inspecting the convergence of the MCMC). Mid- 
dle and Bottom, Panels: Behaviour of model parameter values as 
a function of the iteration number. 




Figure 3. Binned marginal posterior PDFs of the model param- 
eters. These are built after assembling the corresponding samples 
(beyond burn-in) depicted in the middle and bottom panels of 
Fig. [2] Based on these marginal posteriors we may then evalu- 
ate the several Bayesian summary statistics for each of the model 
parameters. 

4.2 On retrieving Do and (r) 

Figs. [4] and [5] result from the application of the ACPS- 
modelling procedure to the full set of realisations of the 
power spectrum (250 realisations for each value of V), al- 
lowing one to assess the potential biases of the Bayesian 
summary statistics of the fitted parameters Do and (F) , re- 
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spectively. The precision associated with these same estima- 
tors is also indicated in the plots. 

When choosing a summary statistic to be used, we 
should make sure that it is representative of the marginal 
probability distribution of the parameter in question, a point 
that all the summary statistics considered seem to satisfy. 
The marginal posterior PDFs of Do a nd (r) are asymmet- 
ric (s ee Fig. [3} and well modelled by (|Fierrv Fraillon et al.l 
1 19981 ): 

f(x) = K(x/sre X p[(x/sr +1 ] , (20) 

where AT is a normalisation constant, s is an adjustable pa- 
rameter describing the shape of the distribution and a de- 
fines the type of distribution (a = for a Boltzmann law and 
high values of a for a Gaussian law) . Furthermore, the set of 
summary parameter values should provide a good fit to the 
data, the quality of which might be assessed by computing 
the rms of the residuals. This latter criterion led us to adopt 
the MAP as our choice of summary statistic. 

It turns out from Figs. [4] and [5] that the different esti- 
mators, although equally robust, present biases whose mag- 
nitude and sign vary with the estimator being considered, a 
direct consequence of the asymmetry of the marginal pos- 
terior PDFs of Do and (r). A way of overcoming these in- 
herent biases would be to present these values normalised to 
the homologous values obtained after performing a similar 
analysis on Sun-as-a-star data covering an equivalent range 
in frequency (i.e. scaling by the ratio of the acoustic cut-off 
frequencies of the Sun and the target in question). 

The MAP summary statistic allows one to determine -Do 
with a relative error of 8 per cent and 17 per cent for V — 8 
and V = 9, respectively. On the other hand, determination 
of (r) is accomplished with a relative error of 13 per cent 
and 18 per cent for V=8 and V = 9, respectively. 

A clear degradation of the precision associated with the 
estimators is seen when considering the set of fainter objects 
(V = 9). This is due to the fact that V is directly coupled to 
the photon shot noise level and hence to S/N. 

Notice that no results have been plotted concerning the 
set characterised by V — 10 as no sensible output was gener- 
ated by the ACPS-modelling procedure. We should mention 
that at this S/N, the automated pipeline is no longer capa- 
ble of supplying an input value for dAiv/dn and hence the 
model p-mode spectrum generated by the ACPS-modelling 
procedure assumes a constant large separation with n. 



4.3 On retrieving i and v s 

With the current experimental setup, we face a scenario 
where the linewidths of the individual modes present in 
the simulated power spectra are systematically larger than 
the rotational splitting, i.e. v s < Y. Consequently, multi- 
plet components are blended together, which strongly cor- 
relates the inclination with the rotational splitting, making 
a successful retrieval of these global parameters not feasi- 
ble. As a further matter, the spectral resolution is too low 
(5v — 0.39 n~PLz) and will only marginally allow the rota- 
tional splitting to be resolved. By again looking at the PDF 
of the inclination for a single realisation of the power spec- 
trum in Fig. [31 we realise that we are in fact retrieving the 
imposed prior since there is no enough evidence in the data 
to proceed otherwise. Fig. [6] corroborates this, if we bear 
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Figure 4. PDFs of the actual errors in the determination of Do, 
which in turn have been normalised after division by their stan- 
dard deviation, tTorr. Several Bayesian summary statistics have 
been considered: a) the MAP, b) the median, c) the mean and d) 
the mode. The solid black line corresponds to the set with V = 8, 
whilst the solid grey line corresponds to the set with V = 9. The 
dashed line represents the standard normal distribution. 
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Figure 5. The same as in Fig.[4]but now regarding the estimation 

of <r>. 
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Figure 6. The posterior mean of the stellar inclination angle, 
i moan , versus the input inclination, i. This plot refers to the set 
of simulations with V = 8. 



in mind that i moan 
r90 



/ O o ip(i\D,I)<ti~ / 0O ip(i\I)di = 
Jqo i sin(i) di — 57° . Similarly, by looking at the PDF of 
the rotational splitting for a single realisation of the power 
spectrum in Fig. [3] we notice that we are in fact retrieving 
the imposed uniform prior, although slightly biased toward 
low values of the splitting. 

In order to demonstrate the full capability of this tool, 
we have applied it to a series of higher S/N realisations of 
the PDS of Boris (assumed to have V = 8), which now cor- 
respond to 6-month long time series (<5!/ = 0.06 /aHz). Since 
we wish to avoid a scenario where v s < V and enter a more 
favourable regime where instead Y < v B < 5vo2, we (i) in- 
creased the reference value for the mean rotational splitting, 
^s Cf £ {2 i/ S Q, 4i/ s q, 6 v s @}, while (ii) simultaneously adopt- 
ing a narrower frequency interval of interest, ranging from 
2100 to 3300 zxHz, at the low-frequency end of the acoustic 
spectrum. This latter step will effectively decrease the refer- 
ence value for the mean linewidth, (r) rcf , since we use solar 
mode linewidths when generating the artificial power density 
spectra, which are known to increase toward higher frequen- 
cies. Finally, we adopted i rcf = 60 ° and imposed uniform 
priors on all the model parameters. Fig. [7] illustrates a case 
where the inclination and the rotational splitting have been 
successfully retrieved, after applying the ACPS-modelling 
procedure to a single realisation of the PDS of Boris. 



5 SUMMARY AND DISCUSSION 

(i) We have developed a new data analysis tool based on 
modelling and fitting the autocovariance of the acoustic PDS 
of a solar-type oscillator. Its main advantage when compared 
to a canonical peak-bagging procedure, relies on the fact 
that there is no need to carry out mode (angular degree I) 
identification prior to performing the fit. This procedure is 
in principle also suitable for being employed in the case of 
evolved stars displaying solar-like oscillations. 

(ii) The implementation of the ACPS-modelling proce- 
dure has been thoroughly described. Furthermore, its auto- 
mated character makes this procedure appropriate for the 
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Figure 7. PDFs of v B and i resulting from the application of the 
ACPS-modelling procedure to a single realisation of the PDS of 
Boris. Reference values are given by v* ei = 4f a Q = 1.6 [iHz and 
i rct = 60°, and are represented by dashed vertical lines. MAP 
estimates are in turn given by t/MAP = 144+0 26 ^jj z anc [ jMAP _ 

55.9 "^Jj'j , therefore being within uncertainties of the reference 
values. 



analysis of a large number of data sets (e.g. arising from the 
Kepler mission). 

(iii) The ACPS conveys information on the large and 
small frequency separations, mode lifetime and rotation. The 
current version of the ACPS-modelling procedure accepts 
Do, u B , (r) and i as free parameters. The prospective in- 
clusion of Av and dA^/dn as additional free parameters is 
envisaged. 

(iv) The tool has been applied to simulated data mim- 
icking what would have been expected for a solar analog 
observed at high cadence during Kepler's survey phase. We 
assessed the potential biases as well as the precision asso- 
ciated with the several Bayesian summary statistics of the 
fitted parameters Do and (r), having been able to decide 
on a preferred summary statistic, namely, the MAP. We 
addressed a way of overcoming the inherent biases. These 
biases would not in any case underpin the usefulness of this 
procedure, since the preferred summary statistic is given to- 
gether with robust estimates of the uncertainties: a MCMC 
sampler is employed in order to obtain the marginal pos- 
teriors for each of the model parameters, from which we 
can estimate the uncertainties. For instance, using the seek 
routine (Quirion et al., in preparation) in order to estimate 
stellar parameters and basing it on a set of asteroseismic 
parameters returned by this procedure, would lead to the 
estimation of a set of stellar p arameters po ssessing robust 
uncertainties (see discussion in lKaroff et al.ll201(t ). 

(v) No sensible output has been generated by the ACPS- 
modelling procedure for the set with V = 10. We could ar- 
gue, of course, that by increasing the effective observational 
length of the simulated data, the method could successfully 
reach a lower S/N since realisation noise is expected to scale 
as 1/VT. 

(vi) Furthermore, retrieval of both i and v s could not be 
achieved with the current experimental setup and a reason 
for that was given based on the intrinsic characteristics of 
the simulated spectra. An example has however been in- 



cluded whereby the full capability of this tool is demon- 
strated. 
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1: procedure PT Metropolis-Hastings 
2: & ,i = ®o , Ki^no 
3: for t = 0, 1, . . . , n;t — 1 do 
4: for i = l,2,..., np do 

5: Propose a new sample to be drawn from a 

proposal distribution: A ~ N(@ t .i\ Ci) 

6: Compute the Metropolis ratio: 

„ _ pWDJkJ) q(Q ti ,|A) 

p(e M |r>,^ 4 ,/) «(A|© M ) 
7: Sample a uniform random variable: 

Ui ~ Uniform(0, 1) 
8: if r > [7i then 

9: t+ i,i = A 

10: else 

ii: e t+ i,i = e M 

12: end if 

13: end for 

14: Ui ~ Uniform(0, 1) 

15: if l/n S wap > ^2 then 

16: Select random chain: 

i ~ UniformInt(l, np — 1) 
17: j = i + 1 

18: Compute r swap : 

_ p(e t ,j|.P,/3,,i>(6>t,,|g,<3 3 ,-0 

19: C/ 3 ~ Uniform(0, 1) 

20: if ^swap ^ U$ then 

21: Swap parameter states of chains i and j: 

&t,i «■ ©tj 
22: end if 

23: end if 

24: end for 

25: return @ M , Vt , i : ft = 1 
26: end procedure 
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Figure Al. Pseudocode- written version of the parallel tempering 
Metropolis— Hastings algorithm, na is the total number of tem- 
pered chains, n- lt is the number of iterations when running the 
MCMC, q(A\&t,i) = N(®t,i\ Cj) is a multivariate Gaussian dis- 
tribution centred on ®t,i and with diagonal covariancc matrix Cj, 
and riawap is the mean number of iterations between successive 
proposals to swap the parameter states of two adjacent chains. 



APPENDIX A: PARALLEL TEMPERING 
METROPOLIS-HASTINGS 



